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Abstract — Joint diagonalisation (JD) is a technique 
used to estimate an average eigenspace of a set of 
matrices. Whilst it has been used successfully in 
many areas to track the evolution of systems via their 
eigenvectors; its application in network analysis is novel. 
The key focus in this paper is the use of JD on matrices 
of spanning trees of a network. This is especially useful 
in the case of real-world contact networks in which a 
single underlying static graph does not exist. The average 
eigenspace may be used to construct a graph which 
represents the 'average spanning tree' of the network or 
a representation of the most common propagation paths. 
We then examine the distribution of deviations from 
the average and find that this distribution in real-world 
contact networks is multi-modal; thus indicating several 
modes in the underlying network. These modes are 
identified and are found to correspond to particular 
times. Thus JD may be used to decompose the behaviour, 
in time, of contact networks and produce average static 
graphs for each time. This may be viewed as a mixture 
between a dynamic and static graph approach to contact 
network analysis. 

Keywords. Social networks, joint diagonalisation, 
graph analysis, spanning tree, human contact networks 

I. Introduction 

Understanding the dynamic structure of contact net- 
works is critical for designing dynamic routing algo- 
rithms [11], epidemic spreading [16] and message 
passing algorithms [10]. 

Time dependent networks are characterised by time 
dependant paths which are characterised by the order 
in which the paths occur. For example, a path between 
3 nodes A ^ B ^ C does not imply a reverse path 
exists; A ^ B ^ C provides no information about 
how C may communicate with A. However, in many 
applications a static graph is constructed which repre- 
sents typically the proportion of time a link was seen 
between two nodes. These static graphs often lose the 
time information which is critical in contact networks. 
However, at a specific time and from a specific node 
there is a single static network representing the paths 
between the root node and the rest of the network. 



JD is used in this paper to look for commonalities 
in these specific static networks and to provide an 
average static network where appropriate. That is, we 
are looking for modes of operation in contact networks. 

Joint Diagonalisation (JD) is a technique that is used 
to track the changes in eigenspace (i.e. eigenvectors 
and eigenvalues) of a system (see Section |llll for exam- 
ples). Eigenvectors and eigenvalues play an important 
role in static network/graph analysis as they can be 
used to determine the centrality of nodes; communities 
and settling times among other things [1]. However, to 
the best of our knowledge tracking eigenspace evolu- 
tion has not been applied to contact/time dependent 
networks previously. This paper examines the use of 
JD in network analysis. 

II. Related work 

Joint diagonalisation has been used in many appli- 
cations where the evolution of a system can be tracked 
smoothly via its eigenspace. For example, Macagnano 
et al. [12] present an algorithm for localisation of 
multiple objects given partial location information. 
As time evolves the location of the objects changes 
smoothly which may be seen through the evolution of 
the eigenvectors of a distance matrix. Other examples 
include blind beam forming [4] and blind source sep- 
aration [20]. 

Sun et al. [19] use tensor analysis to examine time 
dependent networks. A tensor is multi-dimensional 
matrix (for example a set of adjacency matrices) which 
are essentially reduced using PCA to a core tensor. 
This technique is similar in spirit to that presented 
here, the difference is that we are looking to reduce a 
set of spanning trees representing propagation through 
a network; propagation information being preserved. 

Scellato et al. [18] examine the different character- 
istics of contact networks as they evolve over time. 
However, the analysis there is based on forming static 
graphs by amalgamating all links seen in an interval 
of time. This may introduce connections which in 
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Fig. 1. A simple graph and its 6 spanning trees. (The numbers represent the root nodes and probability of observing the tree ex: 1,6,2,3 
are the root nodes for the third tree and this tree is observed with probability 4/14) 



fact are unordered. Graph measures (e.g. the clustering 
coefficient) are then measured from these graphs and a 
time series analysis of these follows. In contrast, here 
the amalgamation is dependent of the contact network 
itself. 

Riolo et al. [16] investigate time dependent epidemic 
networks with a view to constructing transmission 
graphs, directed graphs which indicate the direction 
of transmission of a disease through a network. While 
the aim in this paper is similar, the methodology used 
is significantly different as they examine one time 
infections in real networks. 

III. Theoretical Background 

We begin by defining snowball sampling which 
consists of selecting a root node randomly in the net- 
work with uniform probability and performing a Breath 
First Search (BFS) from this node (i.e. determining a 
set of shortest paths from the source node to every 
other node in the network). This produces a spanning 
tree, H, where is a subset of the original graph 
G{V^E), where V and E denote the vertex and edge 
sets respectively, and \V\ = N denotes the number of 
nodes. We call the starting node the observer or root 
and H, the sample. Figure [T] shows a simple graph 
which will be used for demonstration purposes. In the 
first sample node 5 is selected at random and a shortest 
path first search results in the first tree in Figure [T] In 
this simple graph there are 6 spanning trees shown in 
Figure [T] Note that the distribution of spanning trees 
in this network are not uniform but biased. That is, 
traffic generated uniformly from each node results in 
non-uniform percolation across the network. 

Next we develop a centrality measure which is 
based on standard eigenvector centrality. Eigenvector 
centrality [14] is defined by letting the centrality of 
node i equal the average of the centrality of all nodes 
connected to it: 
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where Aij is element i^j of the (possibly weighted) 
adjacency matrix and A is the largest eigenvalue of 
Ai_^ as can be seen by rewriting Equation [T] in matrix 
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The eigenvector corresponding to 
eigenvector centrality of node i. 

Given M samples of a network, Hi...Hm, the 
question now arises; how can these be combined to 
give a matrix that reflects the sampling bias. We 
propose using a method known as joint diagonalisation 
which produces an average eigenspace of the samples. 
Specifically, we seek an orthogonal matrix such that: 



Hi = UCiU^ \Ji 



(3) 



If U corresponds to the eigenvectors of Hi then Ci 
is diagonal however no matrix U exists in which all 
Ci are diagonal (except for the trivial case in which 
all Hi are equal). Joint diagonalisation seeks average 
eigenvectors U such that the sum of squares of the off 
diagonal elements of Ci are minimised. Specifically: 



M 



U = ar gmin off 2 ( Ci ) 



(4) 



where off 2 is the sum of the off diagonal elements 
squared, called the deviation of Hi from H, di \ 
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where Cf'^ is the k*^ row and j^^ column of Ci. As 
shown in [21] and [3] Equation |4] may be minimised 
efficiently by a sequence of Givens rotations; conver- 
gence and stability properties are proven in [2]. 

Given the average eigenstructure of the sample ma- 
trices an average sampling matrix may be constructed 
from the eigenvector decomposition as: 



H = UCU^ 



(6) 
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Fig. 2. The average sampling graph corresponding to H. 



Where ^ is a matrix in which the entries represent 
the average weight of the Hnks as observed by the 
samples in the network (in a least squares sense) and 
C is the average of diagonals of Hi projected onto U ; 
i.e. the average eigenvalues. A sample based centrality 
may then be constructed from H using the standard 
eigenvector centrality; i.e. by using the eigenvector of 
H corresponding to the maximum eigenvalue. 

A. Simple examples. 

Using the graph shown in Figure [H 100 samples 
are taken by randomly choosing a root node and 
constructing a spanning tree from each. These are 
then jointly diagonalised producing the eigenvectors 
shown in Table HI The sampling centrality is the first 
eigenvector (column 1). 

TABLE I 

The average eigenvector, U from the graph in Figure[T] 
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2 


3 


4 


5 


6 


7 


0.3095 


0.74 


-0.34 


-0.29 


0.13 


-0.35 


0.14 


0.6780 


-0.00 


0.70 


-0.19 


-0.09 


0.00 


-0.09 


0.3370 


0.00 


-0.26 


0.47 


-0.73 


0.00 


0.26 


0.3357 


0.00 


-0.26 


0.45 


0.28 


0.00 


-0.73 


0.3095 


-0.07 


-0.34 


-0.29 


0.13 


0.81 


0.14 


0.1634 


-0.00 


0.14 


0.54 


0.57 


0.00 


0.58 




Fig. 3. The simple graph used in example 2. 




Fig. 4. if for example 2 (Note: the red box contains the node 
number and the sampling centrality). 

It is instructive to view the eigenvector reconstruc- 
tion of H which is formed from the eigenvectors in 
Table U via Equation ^ H is drawn in Figure O Note 
that H is 3. complete weighted graph; weights assigned 
to non-existent links are low and are a consequence of 
taking an average of many graphs. The edge between 
nodes 3 and 4 has a weight of 0.7857 ^ j^, i.e. the 
proportion of trees that use that link (see Figure [T]). 
H represents the sample biased weight of this link; a 
nice result. There are two reasons why these numbers 
are not exactly the same; the first is that, as always, 
there is a slight error introduced when using empirical 
sampling. The second is that the H is based on the 
average eigenspace of a set of sampling trees; this is 
not the same as simply taking the proportion of times 
a link has been observed. 




Fig. 5. if for example 2 with a preference for routes 6^4. 

The second example deals with a preferred route. 
Figure [3] shows a graph with a bridge formed from 



3 



nodes 6<->4 and 5f->3 ; these nodes are critical in 
joining the two parts of the graph. 

Figure |4] shows the weighted graph created by sam- 
pHng as aboveQ. The second stage involves creating 
a preferred route by removing 80% of trees that use 
the link between 5^3. The resultant average graph is 
shown in Figure [51 As can be seen the weight attached 
to link 5^3 is greatly reduced (from 0.9 to 0.2). 

Thus far we have only dealt with samples taken 
on static graphs. However, joint diagonalisation is 
particularly well suited to contact networks. In these 
networks there is no underlying static graph as such, 
but rather a set of contacts that are time dependent. By 
flooding these networks spanning trees may be formed 
and combined by the use of JD. The next section details 
these real- world data sets. 

IV. Data set details 

In this paper, we use four experimental datasets 
gathered by the Haggle Project [7], referred to as Cam- 
bridge, Infocom06; one dataset from the MIT Reality 
Mining Project [6], referred to as MIT. Previously, the 
characteristics of these datasets such as inter-contact 
and contact distribution have been explored in several 
studies [9], to which we refer the reader for further 
background information. These three datasets cover a 
rich diversity of environments, ranging from a quiet 
university town (Cambridge), with an experimental 
period from a few days (Infocom06) to one month 
(MIT). 



Experimental data set 


Cambridge 


Infocom06 


MIT 


Device 


iMote 


iMote 


Phone 


Network type 


Bluetooth 


Bluetooth 


Bluetooth 


Duration (days) 


11 


3 


246 


Granularity (seconds) 


600 


120 


300 


Number of Devices 


36 


78 


97 


Number of contacts 


10,873 


191,336 


54,667 


Average # Contacts/pair/day 


0.345 


6.7 


0.024 



TABLE II 

Characteristics of experimental data sets 



by doing Bluetooth device discovery every five 
minutes. 1 month of data is used here to maintain 
the consistency of users. 
The three experiments are summarised in Table HIl 

and the trace data can be downloaded at CRAWDAD 

database [5]. 

V. Results 

The results below examine the 3 data sets separately 
highlighting the features of each. Finally a synthetic 
contact network with known characteristics is con- 
structed and JD used to extract these characteristics. 




Fig. 6. A typical flooding tree; as seen from node 20. 

Next the distribution of the sample start times is 
examined Figure [51 As can be seen the 5 modes 
correspond to different times in the data set. Modes 
1 and 5 cover the first half of the data while modes 2 
then 3 and then 4 become dominant in that succession. 

A. Cambridge data 

Figure [6] shows a typical sampling tree for the 
Cambridge data set. Node 20 initiates a message and 
it is passed around the contact network; first to nodes 



• In Cambridge, the iMotes were distributed mainly 
to two groups of students from University of 
Cambridge Computer Laboratory, specifically un- 
dergraduate yearl and year2 students, and also 
some PhD and Masters students. This dataset 
covers 11 days. 

• In Infocom06, the trace contains 78 participants. 
Among 78 participants, 34 form 4 subgroups by 
academic affiliations. 

• In MIT, 100 smart phones were deployed to stu- 
dents and staff at MIT over a period of 9 months. 
These phones were running software that logged 
contacts with other Bluetooth enabled devices 

^ Links with low weights are removed for clarity. 
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Community based on Fiedler clustering (Cambridge data 
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Fig. 8. Distribution of 5i (Cambridge data set; kernel smoothing 
is employed for the overall average.). 



Fig. 11. Mean number of nodes susceptible to disease after time t 
(Root infection starts at time 250. Inset focuses on the start of the 
infection.). 




Fig. 9. Distribution of times by mode. 



16 and 6 and from there to the rest of the network. For 
this experiment ten thousand such trees are generatecH 
with the messages starting at a random times and from 
a random node (uniformly distributed). These are then 
combined using Joint DiagonaHzation to form H. 

The average graph, for this data set is represented 
in Figure fTOt a). This representation shows all links in 
the weighted shortest paths of As can be seen the 
nodes split into two groups as expected. These groups 
may be represented by a standard dendrogram based on 
Fiedler vector clustering [8] as shown in Figure |7l The 
groups seen here correspond closely to those found on 
the same data set in [23]. 

The results so far have examined the average be- 
haviour of the contact network which is interesting in 

large sample size is used here to negate random effects. 
However, similar results are found for much smaller sample sizes. 

^We found this to be the clearest means of representing a complete 
weighted graphs. 



itself. However, by examining the distribution of devia- 
tions, 5i, from the average a more interesting behaviour 
may be observed. Figure [8] shows the distribution of 
(^i, z = 1 . . . 10, 000. As can be seen the distribution is 
multi-modal; i.e. the underlying process/contact net- 
work has different modes of operation. A Gaussian 
mixture model [13] is used to determine the different 
modes as shown in Figure [8l 5 different modes are 
identified. This is particularly useful as it allows 
the network to be characterized by different modes of 
behaviour at different times. H for the overall data set 
and for each mode are shown in Figure [TOl Mode one 
shows a highly structured network corresponding to the 
day when the groups are well defined by class year (i.e. 
year 1 and year 2). This structure then becomes less 
well defined as time moves on. Mode 5 is particularly 
interesting as there is an obvious bridge formed by 
nodes 3 and 20. This mode covers the night time and 
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Fig. 12. Graph of shortest paths in H; overall. (MIT). 
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(a) Overall (h) Mode 1 \ ' ^ (c) Mode 2 




Fig. 10. Graph of shortest paths in H for overall and 5 modes. (Cambridge data set; the size of a node is proportional to the sum of 
weights incident on that node) 
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Fig. 13. Community based on Fiedler clustering (MIT data set). 

nodes 3 and 20 are possibly staff who interact with the 
students in the morning. Using H as an indicator, this 
impHes that a disease spread at this time from nodes 
3 and 20 should have the fastest infection rate. Note 
that mode 1 is still dominant in this period; this mode 
is essentially being suspended overnight (due to few 
contacts) with the spanning trees being completed in 
the morning. 

To test the infection rate, an SIR model is con- 
structecQ and a disease is spread through the contact 
network starting at time index 250. The simulation 
is repeated 30 times for each node and the results 
bootstrapped to give estimates of the mean number of 
people susceptible (i.e. those that have not received the 
disease) at time, t, S{t). Figure [TT] shows the results 
of these simulations and as can be seen the number 

"^We cannot be sure as the data has been anonymised. 
^Probability of infection 0.5; infection time Poisson distributed 
with mean 80 time steps, 800 mins. 




Fig. 14. Distribution of 6i (MIT). 




Fig. 15. Distribution of times by mode. (MIT) 

of susceptible people falls most rapidly for infections 
started at nodes 3 and 20, as expected. 

B. MIT data 

The results from the MIT data set show a different 
type of behaviour. There are two main groups (Fig- 
ure [121) the largest of which can be further subdivided 
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Fig. 16. Infocom06 results, (a) Distribution of 6i. (b) Distribution of times by mode, (c-f) H for each mode. 



into three smaller groups (Figure fT3l). 

The distribution of Si is shown in Figure [141 and 
two main modes are identified from this. The MIT 
data set spans a month of data and recurring patterns 
emerge from the data as shown in Figure \T5\ This 
is particularly interesting as it introduces the concept 
of being able to forecast the behaviour of a group at 
regular intervals and design strategies for those specific 
modes. 

C. INFOCOM V6 data 

The Infocom data is summarised in Figure [161 and 
follows the behaviour typically expected at a confer- 



ence. Four modes are identified (Figure [ISa)) which 
correspond to four periods in time (Figure [T6l b)). The 
first mode to occur is mode 3 showing much mixing 
between the delegates (Figure [T6te)). This is probably 
the delegates meeting for coffee before the conference 
begins. This is then followed by two periods of struc- 
tured graphs (i.e. presentations; Figures[T6tc,d)) ending 
with a period of mixing (Figures [T6l f)). 

D. Synthetic contact network 

The first network created in this section is a purely 
random contact network in which 5% of 50 nodes 
are connected at random in each time step. Figure [TSl 
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Fig. 17. (a) The generating graph for submode 2*. (b) H shortest path graph for submode 2*. (c) A spy plot of the generating graph (.) 
and H (o)**. (Synthetic data) * The size of a node is proportional to the sum of weights incident on that node. ** The rows 
and columns have been permuted using approximate minimum degree to highlight the preferential attachment community 
structure in the graph. H has been thresholded using a value of 0.1 < 0.1 i — > 0). 



shows the distribution of Si for this network is uni- 
modal as expected; there is only one underlying pro- 
cess. The distribution also follows a F distributiorEI. 

The second network is more complicated and in- 
volves generating four different behaviours for a con- 
tact network termed generators. A generator consists of 
a static underlying topology representing a set of possi- 
ble contacts. These links are transformed into contacts 
by using a Levy walk (as justified in [9] [15]); a set 
of times are generated from a power law distribution 
and used to demarcate when a contact takes place. 
The generator used is switched every 700 time units as 
shown by the mode indicator in Figure [20l Specifically, 
the first generator employs a Waxman topology [22] 
(a = 0.5,/3 = O.sjZl. The second generator is also a 

^Si is 3. squared quantity which should follow a similar distribution 
to a sample variance; i.e. 6i ~ x^. The F distribution is a 
generalisation of the distribution and so is used. 

^Waxman topology: p{u ^ v) = ae~^^ where a and ^ are 
parameters of the model. The nodes are distributed randomly on a 
grid and the distance between them is d. 



Waxman model with a = 0.7,/3 = 0.3. The third and 
fourth generators are (GLP) generalised linear prefer- 
ential topologies based on preferential attachment [17]. 

As can be seen 3 modes are detected in the data 
(Figure [19]). These correspond with the generator times 
for 2 of the modes (Figure [20]). However, mode 3 
incorporates both generator 3 and generator 4. This 
occurs as generator 3 and 4 are quite similar (both 
based on GLP). 

The samples in mode 3 may be examined separately 
using JD to produce the submodes seen in Figure [21] 
occurring at the times seen in Figure [21] As can be 
seen these submodes are generators 3 and 4. Thus the 
algorithm has successfully recovered the modes in the 
data. At this point we make a note on the transition 
between the modes. It is interesting that this transition 
is not crisp even though the switching between modes 
is. This is because a spanning tree may begin in one 
mode but the message may end in the next mode. 
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Fig. 18. Distribution of Si (Random network). 



Fig. 19. Distribution of 5i (Synthetic). 
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Fig. 20. Distribution of times by mode. (Synthetic) 
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Fig. 21. Distribution of 5i\ submode. (Synthetic). 

VI. Conclusions 

This paper presented a method for extracting dif- 
ferent modes of operation for contact networks. In the 
real- world contact networks examined, several interest- 
ing features where extracted including detection of a 
bridge in the Cambridge data set. The MIT data set in 
contrast, showed a repetitive behaviour which is useful 
for prediction of network behaviour; for example in 
advance of an infection. The INFOCOM data set 
clearly showed the behaviour typical of a conference. 
In producing an average graph based on samples of a 
network, the order of contacts has been preserved and 
in addition the correlation between contacts has been 
preserved. For example aggregation based purely on 
counting the number of times a link is present does 
not take into account the fact that links may typically 
be present together, i.e. the time based correlation be- 
tween links. By using spanning trees the methodology 
takes advantage of a sampling mechanism present in 
many real- world networks; it might not be possible to 
record all contacts but it is often possible to flood a 
message in a network and record the paths taken. It 
is hoped that in future this technique will aid in the 
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Fig. 22. Distribution of times by mode; submode. (Synthetic) 

design of time specific algorithms for time dependent 
networks. 
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